%

%step 0: parameter initialization 
Tinter = 10; % the time interval between frames(ms)
Dthresh = 160; % the threshold of distance between frame in nm
R = 3000;         % the maxium frame allowed to close the gap
Zf = 0; % weight of Z-coordinates
If = 0; % weight of intensity
Zcali = 0.75; % mismatch in z

% step 1: prepare the coordinates structure from 3D-thunderSTORM result
% The format of the result is assumed to be:
%  'id','frame','x [nm]','y [nm]','z [nm]','sigma1 [nm]','sigma2 [nm]','intensity [photon]','offset [photon]','bkgstd [photon]','chi2','uncertainty_xy [nm]','uncertainty_z [nm]'
% setup with the parameters needed
[filename pathname] = uigetfile('.csv', 'select result from thunderstorm processing','multiselect','on');
for ii = 1:length(filename)
    filecurr = filename{ii};
    fileroot = filecurr(1:find(filecurr=='.')-1);
    filesave = ['Coord-' fileroot '.mat'];
    
    fid = fopen([pathname filename{ii}]);
    delimiter = ',';
    startRow = 2;
    formatSpec = '%f%f%f%f%f%f%f%f%f%[^\n\r]';
    ascii_data = textscan(fid, formatSpec, 'Delimiter', delimiter, 'EmptyValue' ,NaN,'HeaderLines' ,startRow-1, 'ReturnOnError', false);
    fclose(fid);
    data = [ascii_data{1:end-1}]; % read the .csv file 
    
    endRow = 1;
    formatSpec = '%q%q%q%q%q%q%q%q%q%q%[^\n\r]';
    fileID = fopen([filename{ii}]);
    dataArray = textscan(fileID, formatSpec, endRow, 'Delimiter', delimiter, 'TextType', 'string', 'ReturnOnError', false, 'EndOfLine', '\r\n');
    txt = [dataArray{1:end-1}];
    fclose(fileID);

    
    Frame = data(:,2); % frame number
    Coordxyz = data(:,3:4); % xy-coordinates
    if strcmp(txt{1,5}, 'z [nm]')
        Coordxyz(:,3) = data(:,5)*Zcali; % z-coordinates
    else
        Coordxyz(:,3) = 0; % no z-coordinates
    end
    Ic = strfind(txt, 'intensity [photon]');
    IndexC = find(not(cellfun('isempty', Ic)));
    Intensity = data(:,IndexC); % spots intensity
    Tn = max(Frame);
    
    for kk = 1:Tn
        Findex =  find(Frame == kk);
        FrameT = Frame(Findex,:);
        CoordT = Coordxyz(Findex,:);
        IntensityT = Intensity(Findex,:);
        id = [1:length(Findex)]';
        Coord = [id,CoordT,IntensityT];
        Coord(:,6:9) = 0;
        Spots(kk).Coord = Coord;
        clear Coord
    end
    
    %step 2: optimize and link the spots
    SpotsLink  = gmtOptimizeBlink( Spots, R, Dthresh, Zf, If ); 
    %step 3: construct the trajectory list
    TrajList  = closeLink(SpotsLink);
    %step 4: construct the tracksFinal sttucture
    tracksFinal = trackConv( TrajList );
    % step 5: save the structures
    save([pathname filesave],'Spots','tracksFinal');
    clear Spots SpotsLink TrajList tracksFinal
end
% clear
% clc
%% analysis 1
% combine the structure from same image
clear
clc
kk = 1;
[filename pathname] = uigetfile('.mat', 'select ttajectories','multiselect','on');
for ii = 1:length(filename)
load([pathname filename{ii}])
 for jj =  1: length(tracksFinal)
     Trace = tracksFinal(jj).Coord;
     if size(Trace,1) > 5
        trackFinal(kk).Coord = Trace;
        kk = kk + 1;
     end
 end
end
tracksFinal = trackFinal;
[filename1 pathname1] = uiputfile('.mat',['Save the trace file of' filename{1}]);
save([pathname1 filename1],'tracksFinal');
%% analysis 2
% gap distribution
gap = [];
[filename pathname] = uigetfile('.mat', 'select ttajectories','multiselect','on');
for ii = 1:length(filename)
load([pathname filename{ii}])
 for jj =  1: length(tracksFinal)
     Trace = tracksFinal(jj).Coord;
     if size(Trace,1) > 1
     G = diff(Trace(:,1));
     gap = [gap;G];
     end
 end
end
hist(gap,[1:max(gap)]);
%% plot the trajectory one by one 3
clear
clc
kk = 1;
[filename pathname] = uigetfile('.mat', 'Select the tracks','multiselect','on');
for ii = 1:length(filename)
 load([pathname filename{ii}])
 filenameRoot = filename{ii};
 TrackROI = tracksRefine.TracksROI;
 for jj =  1: length(TrackROI)
     Trace = TrackROI(jj).Coordinates;
     filenameP = ['3D-' filenameRoot(1:find(filenameRoot=='.')-1) '-' num2str(jj) '.jpg'];
     if size(Trace,1) > 20
         kk = kk + 1
     h =  subplot(2,2,1)
         plot(Trace(:,2)*0.16,Trace(:,3)*0.16,'.-','LineWidth',1);
         title(['X-Y-' filenameRoot num2str(jj)]);
         axis equal
         subplot(2,2,2)
         plot(Trace(:,2)*0.16,Trace(:,4)*0.16,'.-','LineWidth',1);
         title('X-Z');
         axis equal
         subplot(2,2,3)
         plot(Trace(:,3)*0.16,Trace(:,4)*0.16,'.-','LineWidth',1);
         title('Y-Z');
         axis equal
     saveas(h,filenameP,'jpg');
     end
 end
end

%% plot the trajectory one by one on the cell 4
clear
clc

[filename1 pathname1] = uigetfile('.tif', 'Select the BFimage');
[filename pathname] = uigetfile('.mat', 'Select the tracks');

 load([pathname filename])
 
ImageBF = imread([pathname1 filename1]);

 TrackROI = tracksRefine.TracksROI;
 h = figure('position',[100,100,1000,1000]);
 imshow(ImageBF,[],'InitialMagnification',200);
 set(h,'position',[100,100,500,500]);
 hold on
 for jj =  1: length(TrackROI)
     Trace = TrackROI(jj).Coordinates;
%      filenameP = ['Direct-' filenameRoot(1:find(filenameRoot=='.')-1) '-' num2str(jj) '.fig'];
%      if size(Trace,1) > 20
     plot(Trace(:,2),Trace(:,3),'.-','LineWidth',1);
         h=subplot(2,2,1)
         plot(Trace(:,2),Trace(:,3),'.-','LineWidth',1);
         title(['X-Y-' num2str(jj)]);
         axis equal
         subplot(2,2,2)
         plot(Trace(:,2),Trace(:,4),'.-','LineWidth',1);
         title('X-Z');
         axis equal
         subplot(2,2,3)
         plot(Trace(:,3),Trace(:,4),'.-','LineWidth',1);
         title('Y-Z');
         axis equal
         waitforbuttonpress
         close all
     title(filenameP);
   saveas(h,filenameP,'fig');
     end
%  end
%% plot the trajectory one by one on the cell 5
clear
clc

[filename1 pathname1] = uigetfile('.tif', 'Select the BFimage');
[filename pathname] = uigetfile('.mat', 'Select the tracks');

 load([pathname filename])
 
ImageBF = imread([pathname1 filename1]);

 TrackROI = tracksRefine.TracksROI;
 h = figure('position',[100,100,1000,1000]);
 imshow(ImageBF,[],'InitialMagnification',200);
 set(h,'position',[100,100,500,500]);
 hold on
 for jj =  1: length(TrackROI)
     Trace = TrackROI(jj).Coordinates;
    filenameP = ['Direct-' filenameRoot(1:find(filenameRoot=='.')-1) '-' num2str(jj) '.fig'];
      if size(Trace,1) > 5
     plot(Trace(:,2),Trace(:,3),'.-','LineWidth',1);
         h=subplot(2,2,1)
         plot(Trace(:,2),Trace(:,3),'.-','LineWidth',1);
         title(['X-Y-' num2str(jj)]);
         axis equal
         subplot(2,2,2)
         plot(Trace(:,2),Trace(:,4),'.-','LineWidth',1);
         title('X-Z');
         axis equal
         subplot(2,2,3)
         plot(Trace(:,3),Trace(:,4),'.-','LineWidth',1);
         title('Y-Z');
         axis equal
         waitforbuttonpress
         close all
     title(filenameP);
   saveas(h,filenameP,'fig');
     end
 end
 
 %% make structure with specific trajectories 6
 clear
 clc
 [filename pathname] = uigetfile('.xlsx','Select the excel file with trace number');
 [data txt] = xlsread([pathname filename]);
 for ii = 1 : size(data,1)
     ROI = data(ii,1);
     number = data(ii,2);
     filename = ['TraceROI' num2str(ROI) '-ref-all.mat'];
     load(filename);
     Track = tracksRefine.TracksROI;
     Struc.TracksROI(ii).Coordinates = Track(number).Coordinates;
 end
 
 